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*H : 1 Introduction 



The Computational Singular Perturbation (CSP) method is one of several 
so-called reduction methods developed in chemistry to systematically decrease 
the size and complexity of systems of chemical kinetics equations. The method 
was first proposed by Lam and Goussis [HI El 13 El E] and is widely used, for 
example, in combustion modeling jH [TUl HU H21 UH UH] • 

The CSP method is generally applicable to systems of nonlinear ordi- 
nary differential equations (ODEs) with simultaneous fast and slow dynamics 
where the long-term dynamics evolve on a low- dimensional slow manifold in 
the phase space. The method is essentially an algorithm to find successive 
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approximations to the slow manifold and match the initial conditions to the 
dynamics on the slow manifold. 

In a previous paper [20], we focused on the slow manifold and the ac- 
curacy of the CSP approximation for fast-slow systems of ODEs. In such 
systems, the ratio of the characteristic fast and slow times is made explicit by 
a small parameter e, and the quality of the approximation can be measured in 
terms of e. By comparing the CSP manifold with the slow manifold found in 
Fenichel's geometric singular perturbation theory (GSPT, [2112]), we showed 
that each application of the CSP algorithm improves the asymptotic accuracy 
of the CSP manifold by one order of e. 

In this paper, we complete the analysis of the CSP method by focusing 
on the fast dynamics. According to Fenichel's theory, the fast-slow systems 
we consider have, besides a slow manifold, a family of fast stable fibers along 
which initial conditions tend toward the slow manifold. The base points of 
these fibers lie on the slow manifold, and the dynamics near the slow manifold 
can be decomposed into a fast contracting component along the fast fibers 
and a slow component governed by the motion of the base points on the slow 
manifold. By comparing the CSP fibers with the tangent spaces of the fast 
fibers at their base points, we show that each application of the CSP algorithm 
also improves the asymptotic accuracy of the CSP fibers by one order of e. 

Summarizing the results of [20] and the present investigation, we conclude 
that the CSP method provides for the simultaneous approximation of the slow 
manifold and the tangents to the fast fibers at their base points. If one is 
interested only in the slow manifold, then it suffices to implement a reduced 
(one-step) version of the algorithm. On the other hand, if one is interested in 
both the slow and fast dynamics, then it is necessary to use the full (two-step) 
CSP algorithm. Moreover, only the full CSP algorithm allows for a linear 
matching of any initial data with the dynamics on the slow manifold. 

This paper is organized as follows. In Section |2] we recall the relevant 
results from Fenichel's theory and set the framework for the CSP method. In 
Section we outline the CSP algorithm and state the main results: Theo- 
rem concerning the approximation of the slow manifold, which is a verbatim 
restatement of |2Til Theorem 3.1]; and Theorem 13 . 21 concerning the approxima- 
tion of the tangent spaces of the fast fibers. The proof of Theorem 13.21 is given 
in Section |3J In Section we revisit the Michaelis-Menten-Henri mechanism 
of enzyme kinetics to illustrate the CSP method and the results of this article. 
Section |U] is devoted to a discussion of methods for linearly projecting initial 
conditions on the slow manifold. 
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2 Slow Manifolds and Fast Fibers 



Consider a general system of ODEs, 

for a vector-valued function x = x(t) e j n a sm0 oth vector field g. 

For the present analysis, we assume that n components of x evolve on a time 
scale characterized by the "fast" time t, while the remaining m components 
evolve on a time scale characterized by the "slow" time r = et, where e is 
a small parameter. (The explicit identification of a small parameter e is not 
necessary for the applicability of the CSP method; a separation of time scales 
is sufficient.) We collect the slow variables in y E R m and the fast variables 
in z G R n . Thus, the system (12. 1|) is equivalent to either the "fast system" 

y' = e 9l (y,z,e), (2.2) 
z' = g 2 (y,z,e), (2.3) 



or the "slow system" 



y = gi(y,z,e), (2.4) 

ez = g 2 {y,z,e). (2.5) 



(A prime ' denotes differentiation with respect to t, a dot ' differentiation 
with respect to r.) The fast system is more appropriate for the short-term 
dynamics, the slow system for the long-term dynamics of the system (J2.1j) . 

In the limit as e tends to 0, the fast system reduces formally to a single 
equation for the fast variable z, 

z' = g 2 (y,z,Q), (2.6) 

where y is a parameter, while the slow system reduces to a differential equation 
for the slow variable y, 

y = gi {y,z,0), (2.7) 
with the algebraic constraint g 2 (y, z, 0) = 0. 

We assume that there exist a compact domain K and a smooth function ho 
defined on K such that 

g 2 (y,h {y),0) = 0, y e K. (2.8) 
The graph of h defines a critical manifold A4o, 

Mo = {{y, z) e R m+n : z = h (y), y e K}, (2.9) 
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and with each point p = (y, ho{y)) E A4o is associated a fast fiber J 7 ^, 

Tl = {(y, z) e R m+n : z e R n }, p e M . (2.10) 

The points of .Mo are fixed points of Eq. (j2.6|l . If the real parts of the eigen- 
values of D z g 2 (y, ho(y), 0) are all negative, as we assume, then .Mo is asymp- 
totically stable, and all solutions on T§ contract exponentially toward p. 

If e is positive but arbitrarily small, Fenichel's theory [21 13] guarantees 
that there exists a function h £ whose graph is a slow manifold A4 £ , 

M e = {(y, z) e R m+n : z = h £ ( y ), y e K}. (2.11) 

This manifold is locally invariant under the system dynamics, and the dynam- 
ics on M. £ are governed by the equation 

y = 9i(y,h £ (y),e), (2.12) 

as long as y £ K. Fenichel's theory also guarantees that there exists an 
invariant family JF £ , 

?e = |J T>, (2.13) 

of fast stable fibers T\ along which solutions relax to Ai £ . The family is invari- 
ant in the sense that, if cj>t denotes the time-t map associated with Eq. (|2.1j) . 
then 

^(ff)cJf w , P eM £ . (2.14) 

The collection of fast fibers J 7 ? foliates a neighborhood of Ai £ . Hence, the 
motion of any point on JF^ decomposes into a fast contracting component 
along the fiber and a slow component governed by the motion of the base 
point of the fiber. Also, M. £ is C(e)-close to .Mo, with 

h e {y) = h {y) + ehi{y) + s 2 h 2 (y) + ■ • ■ , e | 0, (2.15) 

and Tl is (9(£)-close to J 7 ^ in any compact neighborhood of Ai £ . 

Remark 2.1. Typically, the manifold M. £ is not unique; there is a family of 
slow manifolds, all having the same asymptotic expansion ()2.15|) to all orders 
in e but differing by exponentially small amounts ((9(e _C//e ), c > ). 



3 The CSP Method 

The CSP method focuses on the dynamics of the vector field g(x), rather than 
on the dynamics of the vector x itself. 
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Writing a single differential equation like (|2.1|) as a system of equations 
amounts to choosing a basis in the vector space. For example, in Eqs. ([2.2)1 - 
([2. 3 J) , the basis consists of the ordered set of unit vectors in R m+n . The 
coordinates of g relative to this basis are eg\ and g-i- If we collect the basis 
vectors in a matrix in the usual way, then we can express the relation between 
g and its coordinates in the form 

»=(»Oft)- (31) 

Note that the basis chosen for this representation is the same at every point 
of the phase space. The CSP method is based on a generalization of this idea, 
where the basis is allowed to vary from point to point, so it can be tailored to 
the local dynamics near Ai £ . 

Suppose that we choose, instead of a fixed basis, a (point-dependent) 
basis A for R m+n . The relation between the vector field g and the vector / of 
its coordinates relative to this basis is 



9 = Af. (3.2) 

Conversely, 

/ = Bg, (3.3) 

where B is the left inverse of A, BA = I on R m+n . In the convention of the 
CSP method, A is a matrix of column vectors (vectors in R m+n ) and B a 
matrix of row vectors (functionals on R m+n ). 

The CSP method focuses on the dynamics of the vector /. Along a 
trajectory of the system (j2.1jl . / satisfies the ODE 

| = A/, (3.4) 

where A is a linear operator (THl E0j , 

A = B(Dg)A +^§A = B(Dg)A - B^ = B[A, g}. (3.5) 

Here, Dg is the Jacobian of g, dB/dt = (DB)g, dA/dt = (DA)g, and [A,g] is 
the Lie bracket of A (taken column by column) and g. The Lie bracket of any 
two vectors a and g is [a,g] = (Dg)a — (Da)g; see [Tlj . 

It is clear from Eq. ()3.4|) that the dynamics of / are governed by A, so 
the CSP method focuses on the structure of A. 
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Remark 3.1. It is useful to see how A transforms under a change of basis. If 
C is an invertible square matrix representing a coordinate transformation in 
R m+ri , and A = AC and B = C^B, then 

A = B(Dg)A-B^ = C- 1 B(Dg)AC-C- 1 B^l 

An 

= C~ l KC -C- 1 ^. (3.6) 
dt v 1 

Hence, A does not transform as a matrix, unless C is constant. 



3.1 Decompositions 

Our goal is to decompose the vector / into its fast and slow components. Sup- 

( f 1 \ 

pose, therefore, that we have a decomposition of this type, / = I ' 2 , where 



. P . 

f 1 and f 2 are of length n and m, respectively, but not necessarily fast and 
slow everywhere. The decomposition suggests corresponding decompositions 

B 1 



of the matrices A and B, namely A = (Ai, A 2 ) and B = y ^ 2 j , where Ai is 

an (m + n) x n matrix, A 2 an (m + n) x m matrix, B 1 an n x (m + n) matrix, 
and -B 2 an m x (m + n) matrix. Then, f 1 = B 1 g and f 2 = B 2 g. 

The decompositions of A and B lead, in turn, to a decomposition of A, 

A 12 W^i,*?] ^[A 2 ,<?]\ , ^ 

v a 21 a 22 J ~ v 5 2 [^i^] 5 2 [a 2 ,^] ; • {6J} 

The off-diagonal blocks A 12 and A 21 are, in general, not zero, so the equations 
governing the evolution of the coordinates f 1 and f 2 are coupled. Conse- 
quently, f 1 and f 2 cannot be identified with the fast and slow coordinates of 
g globally along trajectories. The objective of the CSP method is to construct 
local coordinate systems (that is, matrices A and B) that lead to a block- 
diagonal structure of A. We will see, in the next section, that such a structure 
is associated with a decomposition in terms of the slow manifold and the fast 
fibers. 

Remark 3.2. Note that the identity BA = I on R m+n implies four identities, 
which are summarized in the matrix identity 

B 1 A 1 B l A 2 \_( In 
B 2 A, B 2 A 2 ) { I m 



(3.8) 
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3.2 Block-Diagonalization of A 



In this section we analyze the properties of A relative to a fast-slow decompo- 
sition of the dynamics near A4 £ . 

Let T P T £ and T P M. £ denote the tangent spaces to the fast fiber and the 
slow manifold, respectively, at the base point p of the fiber on hA £ . (Note 
that <X\mT p T £ = n and (\m\T p M. £ = m.) These two linear spaces intersect 
transversally, because Ai £ is normally hyperbolic and compact, so 

R m +" = 7 P T £ © T P M £ , peM £ . (3.9) 

Let Af be an (m+n) xn matrix whose columns form a basis for T p J : £ and A s an 
(m+n) x m matrix whose columns form a basis for T p Ai £ , and let A = {Af, A s ). 
(We omit the subscript p.) Then A is a (point-dependent) basis for R m+n that 
respects the decomposition (|3.9j) . We recall that 1 M. £ = UpeA^fe %>M £ ) and 
TT £ = UpgA4 e (P; 1pJ~e) are the tangent bundles of the slow manifold and the 
family of the fast fibers, respectively. (A general treatment of tangent bundles 
of manifolds is given in [IJ Section 1.7].) 

The decomposition ()3.9|) induces a dual decomposition, 

R m+n = NpM ^ M ^ p G M ^ (31Q) 

where N P M. £ and M P T £ are the duals of T P M. £ and T P T £ , respectively, in 
R m+n . (Note that d\mj\f p M. £ = n and d\mj\f p T £ = m.) The corresponding 

decomposition of B is B = ( yj_ j , where the rows of i? s_L form a basis for 

M P M. £ and the rows of -B^ -1 a basis for M P T £ . Furthermore, 



B s± A f B S± A S \ f I n 
B^A f B^A S ) { I m 



(3.11) 



The decompositions of A and B lead, in turn, to a decomposition of A, 

/ B^[A f ,g] B^[A s ,g] \ 
A -{B^[A f ,g] Bf±[A.,g)J- ^ 

This decomposition is similar to, but different from, the decomposition (|3.7|) . 
The following lemma shows that its off-diagonal blocks are zero. 



Lemma 3.1 The off- diagonal blocks in the representation of A are zero 

at each point p 6 M. £ . 
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Proof. Since B S± A S = on M. £ and A4 £ is invariant, we have 

j t {B S± A S ) = D(B s± A s )g = (DB sL )(g, A.) + B s± ((DA s )g) = 0. (3.13) 

(DB S± is a symmetric bilinear form; its action on a matrix must be understood 
as column-wise action.) 

Also, g G TM e , so B s± g = on M. e . Hence, the directional derivative 
along A s (taken column by column) at points on M. s also vanishes, 

D(B a± g)A 8 = (DB s± )(A s ,g) + B s± (Dg)A s = 0. (3.14) 

Subtracting Eq. ()3.13|) from Eq. (|3.14j) . we obtain the identity 

B s± [A s ,g] = B s± ((Dg) A s - (DA s )g) = 0. (3.15) 

The proof for the lower left block is more involved, since the fast fibers are 
invariant as a family. Assume that the fiber T$ at p G A4 e is given implicitly 
by the equation F(q; p) = 0, q G T^. Then the rows of (D q F) (q; p) form a basis 
for N q T e , so there exists an invertible matrix C such that B^ = C(D q F). 

Since the rows of (D q F)(q;p) span N q T e , we have (D q F)(q;p)Af(q) = 0. 
This identity holds, in particular, along solutions of (|2.1|) . so 

j t ((D q F)(q;p)A f (q)) = ((D* q F)(q;p)) (g(q), A f (q)) 

+ ((D pq F)(q;p))(g(p),A f (q)) 
+ ((D q F)(q;p))(DA f (q))g(q) 
= 0. (3.16) 

The family of the fast fibers is invariant under the flow associated with (J2.1j) . 
so if F(q;p) = 0, then also F(q(t);p(t)) = and, hence, 

= ((D q F)(q;p))g(q) + ((D p F)(q; p)) g(p) = 0. (3.17) 

Next, we take the directional derivative of both members of this equation along 
Af, keeping in mind that (Dg)(p)Af(q) = because the base point p does not 
vary along Af. (Recall that the columns of Af(q) span T q J r e .) We find 

((D 2 q F)(q;p)) (A f (q),g(q)) + ((D q F)(q;p)) (Dg(q))A f (q) 

+ ((D pq F)(q;p))(A f (q),g(p)) = 0. (3.18) 

But the bilinear forms D q F and D pq F are symmetric, so subtracting Eq. (|3.1fij) 
from Eq. ()3.18|) and letting q = p, we obtain the identity 

(D q F)(p;p) {{Dg)A f - (DA f )g) (p) = 0. (3.19) 
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Hence, B f± [A f ,g}(p) = C{D q F)(p;p)[A f , g}(p) = 0, and the proof of the 
lemma is complete. I 



The lemma implies that the representation (j3.12J) is block-diagonal, 

f B^[A f ,g] \ 

Consequently, the decomposition ()3.9|) reduces A. In summary, if we can 
construct bases Af and A s , then we will have achieved a representation of A 
where the fast and slow components remain separated at all times and the 
designation of fast and slow takes on a global meaning. 



3.3 The CSP Algorithm 

The CSP method is a constructive algorithm to approximate Af and A s . One 
typically initializes the algorithm with a constant matrix A^°\ 

A«» = (a?\ Af) = ( < J \ . (3.21) 

Here, Af± is an m x n matrix, annxm matrix, and the off-diagonal blocks 
A12 and are full-rank square matrices of order m and n, respectively. A 



common choice is A^ = 0. We follow this convention and assume, henceforth, 
that A ( S = 0, 

(o) " 



A»> = (A?\ Af) = ( ^ ) . (3.22) 



(Other choices are discussed in |20j.) The left inverse of A^°> is 
Bn» = ( B J?A = ( 3 3 



(0) " \ Bf 0) ) ~ \ B% 



(3.23) 



The algorithm proceeds iteratively. For q — 0,1, ... , one first defines the 
operator A( 9 ) in accordance with Eq. ()3.5|) . 

A (9) = - = ( $f $g ) , (3.24) 
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and matrices U( q ) and L( g ), 

^-U o J' L w - Us^g))- 1 ° J' ( } 

Then one updates and Bi q \ according to the formulas 

A <s+i) = A^(I-U iq) )(I + L {q) ), (3.26) 
B {q+1) = (I-L {q) )(I + U {q) )B (q) , (3.27) 

and returns to Eq. ()3.24|) for the next iteration. 

Remark 3.3. Lam and Goussis jH] perform the update ()3.26J) - (j3.27J) in two 

steps. The first step corresponds to the postmultiplication of A^ q ' with / — Um 
and premultiplication of B^ with / + U( q ), the second step to the subsequent 
postmultiplication of A^(I — Ut q )) with / + Lr q ) and premultiplication of (/ + 
U( q ))B( q ) with I - 



3.4 Approximation of the Slow Manifold 

After q iterations, the CSP condition 

£ ( > = 0, 9 = 0,1,..., (3.28) 

identifies those points where the fast amplitudes vanish with respect to the 
then current basis. These points define a manifold that is an approximation 
for the slow manifold Ai £ . 

For q = 0, .BL is constant and given by Eq. (j3.23|) . Hence, the CSP 
condition (|3.28|) reduces to the constraint g 2 (y,z,e) = 0. In general, this 
constraint is satisfied by a function z = i/jr \(y,s). The graph of this function 

defines K,f \ the CSP manifold (CSPM) of order zero. Since the constraint 
reduces at leading order to the equation g 2 (y,z, 0) = 0, which is satisfied by 
the function z = h (y), Ke may be chosen to coincide with A4q to leading 
order; see Eq. (|2.9j) . 

For q — 1,2,..., the CSP condition takes the form 

B ( q )(y^( q -i)(y^)^)9(y,z,e) = o, g = i,2, .... (3.29) 

The condition is satisfied by a function z = ip( q )(y,E), and the manifold 

= {(y, z):z = if> (q) (y, e), y G K}, q = 0, 1, . . . (3.30) 
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defines the CSP manifold (CSPM) of order q, which is an approximation of 
Ai £ . The following theorem regarding the quality of the approximation was 
proven in |20j . 

Theorem 3.1 \2(A Theorem 3.1] The asymptotic expansions of the CSP man- 
ifold K.i and the slow manifold Ai £ agree up to and including terms ofO(e q ), 

Q 

iJ iq) (-,e) = J2^h, + 0(e C!+1 ), £ |0, g = 0,l,.... (3.31) 

j=0 

3.5 Approximation of the Fast Fibers 

We now turn our attention to the fast fibers. The columns of Af(y, h e (y)) span 
the tangent space to the fast fiber with base point p = (y, h £ (y)), so we expect 
that Ai defines an approximation for the same space after q applications of 
the CSP algorithm. We denote this approximation by df\y) and refer to it 
as the CSP fiber (CSPF) of order q at p, 

C®(y) = span (cols (A?(y, ^ q) (y, e), e))). (3.32) 

We will shortly estimate the asymptotic accuracy of the approximation, but 
before doing so we need to make an important observation. 

Each application of the CSP algorithm involves two steps, see Remark 3.3. 
The first step involves U and serves to push the order of magnitude of the upper 
right block of A up by one, the second step involves L and serves the same 
purpose for the lower left block. The two steps are consecutive. At the first 
step of the qth iteration, one evaluates Bis on K!f~ 1 ^ to find fC^f 1 by solving 
the CSP condition ()3.28|) for the function ifj^. One then uses this expression 
in the second step to update A and B, thus effectively evaluating A± on K,^ 
rather than on /Ce 9_1 ' ) . 

The following theorem contains our main result. 

Theorem 3.2 The asymptotic expansions of d (y) and T p T e , where p = 
(y, h e (y)) G Ai e , agree up to and including terms of 0(e q ), for all y e K and 
for q = 0,1, .... 

Theorem 13 . 21 implies that the family df 1 = [j p£Ms {p, df\y)) is an 0(e q )- 
approximation to the tangent bundle TT £ . 
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The proof of Theorem 13.21 is given in Section 0] The essential idea is to 
show that, at each iteration, the asymptotic order of the off-diagonal blocks 
of A( 9 ) increases by one and and B?^ become fast and fast -1 , respectively, 
to one higher order. As a consequence, in the limit as q — > oo, A( 9 ) — > A, 
^4(9) _> anc [ 5, ^ _> where A, A, and £> are ideal in the sense described 
in Section EUl 

Remark 3.4. If, in the second step of the CSP algorithm, A^ were evaluated 
on K}§ ^ instead of on K}i\ the approximation of TT e might be only 0(e q ~ 1 )- 
accurate. However, see Section |S] for an example where the approximation is 
still 0(e q ). 



4 Proof of Theorem 13.2 

The proof of Theorem 13.21 is by induction on q. Section |4~T1 contains an auxil- 
iary lemma that shows that each successive application of the CSP algorithm 
pushes A closer to block-diagonal form. The induction hypothesis is formu- 
lated in Section l4~2l the hypothesis is shown to be true for q = in Section IP) 
and the induction step is taken in Section 14.41 



4.1 Asymptotic Estimates of A 

As stated in Section El the goal of the CSP method is to reduce A to block- 
diagonal form. This goal is approached by the repeated application of a two- 
step algorithm. As shown in [20], the first step of the algorithm is engineered 
so that each application increases the asymptotic accuracy of the upper-right 
block AH by one order of e; in particular, A^ = 0(e q ) on /Ci 9) [2H Eq. (5.25)]. 
We now complete the picture and show that each application of the second 
step increases the asymptotic accuracy of the lower-left block A?\ by one order 
of e, when the information obtained in the first step of the same iteration is 
used. In particular, A?\ = 0(e q+1 ) on K.i q+1 \ where JCi q+1 ^ has been obtained 
in the first step of the (q + l)th refinement. 



Lemma 4.1 For q = 0, 1, . . 



when A( g ) is evaluated on K, 



(9+1) 

£ 
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Proof. The proof is by induction. The desired estimates of A}\, A^, and 

AJJ on JCi q) were established in (213 Eqs. (5.24), (5.25), (5.27)]. Since the 

asymptotic expansions of K!f +l ^ and differ only at terms of 0(e q+1 ) or 
higher ([201 Theorem 3.1]), these estimates of A^, AR, and A^ are true also 

on ldf +l \ It only remains to estimate A?\. 

Consider the case q = 0. Let A^ ^ be the coefficient of e J in the asymptotic 

expansion of A^Jy, ip^(y), e). The estimate A?i = O(e) on /Ce follows if we 
can show that A^ ^ = 0. It is already stated in [201 Eq. (4.30)] that A^^ = 

on Ke ■ Furthermore, [213 Theorem 3.1] implies that the asymptotic expan- 
sions of ip(i) and ^(o) agree to leading order. Thus, the asymptotic expansions 
of A?Qv(y, ^)(o)(y), e) and A?A(y, ip^(y), e) also agree to leading order, and the 
result follows. 

Now, assume that the asymptotic estimate holds for 0,1,... , q. From 
Eq. ()3.6|) we obtain 

A S) - L (<?) A (g) + A (J L (<?) - L (q) A }q) L (q) ~ A H) U (q) L (q) 

~ L (<i) U (<i) A tl) + L (q) A H) U (q) L (g) - L (q) U (i) A2 (q) L (i) 
+ kq) U (q) A ti) U (q) L M + ( D h 9 )) 9 + hi) (( DU (i)) 9) L [q) . (4.2) 

The first two terms in the right member sum to zero, by virtue of the defini- 
tion (j3.25|) of The next seven terms are all 0(e q+2 ) or higher, by virtue 
of the induction hypothesis. Finally, the last two terms are also 0(e q+2 ) or 
higher, by the induction hypothesis and [213 Lemma A. 2]. | 



A 21 

A (9+l) ~ 



4.2 The Induction Hypothesis 

The CSPF of order q, Ce(y), is defined in Eq. f|3.32|) to be the linear space 
spanned by the columns of the fast component, A± (y, xp! q \ : e) : of the basis A^ q \ 
Thus, to prove Theorem 13 .21 it suffices to show that the asymptotic expansions 
of A^\y, ip^q), s) and the space tangent to the fast fiber, T v T e , agree up to and 

including terms of 0(e q ), for p = (y, h e (y)) and for q = 0, 1, The central 

idea of the proof is to show that each successive application of the CSP method 
pushes the projection of on TM. £ along TT e to one higher order in e. 

We express A^ q \ generated after q applications of the CSP algorithm, in 
terms of the basis A, 

A( q \y,z,e) = A(y,h £ ,e)Q( q \y,z,e), q = 0, 1, . . . . (4.3) 
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Since Br g \ and B are the left inverses of and A, respectively, we also have 
B( g )(y,z,e) = R( q) (y,z,e)B(y,h e ,e), q = 0, 1, . . . , (4.4) 
where = {Q^)~ l . Introducing the block structure of and R( q ), 

we rewrite Eqs. f!4.3|) and (|4.4|) as 

4 9) = ^?+a,oS, 4* ) = ^/^+aoS? ) (4.6) 

and 

B U = K) ±BS± + <) ±s/± > 4» = 4t BS± + R lt) ±Bf± > ( 4 - 7 ) 

for q = 0, 1, 

Equation ()4.7j) shows that is the projection of A± on TM. £ . Thus, 

to establish Theorem l3.2l we only need to prove the asymptotic estimate Q^j = 
0(e q+1 ). The proof is by induction on q, where the induction hypothesis is 

= (^i) ^(t))' 9 = 0,1,.... (4.9) 

Remark 4.1. Although the estimate of Q^j is sufficient to establish Theo- 
rem EiH we provide the estimates of all the blocks in Eqs. ()4.8j) - (j4.9j) because 
they will be required in the induction step. 

The validity of Eqs. ()4.8j) - (j4.9)) for q = is shown in Section 14.31 The 
induction step is carried out in Section POl 



4.3 Proof of Theorem EH for q = 

We fix q = and verify the induction hypothesis for Q^> and Rm- By Eq. (|4.3jl 

q(°) = BA<®, (4.10) 

whence 

14 



It suffices to show that the lower-left block is zero to leading order, since the 
other blocks are all 0(1). We do this by showing that Q2* 0) = 0. By Eq. I[PT)) . 

QT=Bl x A?' \ (4.12) 



Bq 1 spans N p J-q for every p e /Ce . Also, z is constant on ApJF , so 5^ = 
(B 1 ?- 1 , 0), where i? 1 ^- 1 is a full-rank matrix of size m. Last, 4°< 0) = 4 0) = 

^(o) )j by Eq. (jH22J). Substituting these expressions for and Af'^ 1 

into Eq. (ETT^ . we obtain that Qff ] = 0. 

The induction hypothesis on can be verified either by a similar argu- 
ment, or by recalling that Rm\ = (Q^) -1 , where was shown above to be 
block-triangular to leading order. 



4.4 Proof of Theorem EH for q = 1, 2, . . . 

We assume that the induction hypothesis (|4.8j) - (|4.9j) holds for 0, 1, . . . , q and 
show that it holds for q + 1. The proof proceeds in four steps. In step 1, we 
derive explicit expressions for R( q +i) and Q^ q+l > in terms of R^ q ) and Q^; these 
expressions also involve U^ q ) and L( g ). In step 2, we derive the leading-order 
asymptotics of Uu), and in step 3 the leading-order asymptotics of L^. Then, 
in step 4, we substitute these results into the expressions derived in step 1 to 
complete the induction. 

Step 1. We derive the expressions for Q (9+1) and R{ q+ i). Equations (|Oj) 
and ()4.4j) . together with the update formulas ()3.26|) for and ()3.27|) for 
B {q) , yield 

Q&r+D = gW(J-C/ (g) )(J + L (?) ), (4.13) 
= (I-L {q) )(I + U ig) )R (q) . (4.14) 

Q { $UwL iq) , (4.15) 

(4.16) 

( 4 - 17 ) 

(4.18) 
(4.19) 



In terms of the constituent blocks, we have 
Qif~ ^ = Q$ + Q^fL^ 
Q*2f 1 = Qy ~ QvjUfa)) 

Q 
Q 

and 

pls_L pls-L i rr E?2s_L 

-"(g+l) - U (q) + ty (9)- K (< ? ) > 



,(9+1) 
1/ 


= <#/' 


+ Qf, 


,(9+1) 
2/ 


— Qy 


-Qf, 


,(9+1) 
Is 


= off 


+ 


,(9+1) 
2s 


= «S? 


- Q'H 
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*(£"!) = R \^ + u («) R ll^ ( 4 - 2 °) 

E>2s_L I?2s_l_ t ols-L r TT Z?2s± /yi r>i\ 

it ( g+i) - tt {q) -L {q) K iq) - L {q) U {q) K {q) , (4.21) 

<i) = <) ± -^)< / ) ± -^)^)<) ± - ( 4 -22) 



Step 2. We derive the leading-order asymptotics of the matrix U( q ). 

Recall that U {q) = (Agj) -1 Ag). Moreover, AJJ is strictly 0(1) and AfJ is 
strictly C?(e«) by Lemma Ol Hence, £/ (g) = £/ (<M )£ 9 + 0(£ 9+1 ), with U {m) = 
(A^ ^) _1 Ag (j ^. Therefore, it suffices to derive the leading order asymptotics 
of these blocks of A. 

By definition, A( 5 ) = B^[A^ q \ g}. Therefore, 

^Uk'j] b U ^\ 9] )- (4 ' 23) 

The individual blocks of are obtained by substituting Eqs. ()4.6|) and (14.71) 
into Eq. (|4.23|) . We observe that one-half of all the terms would vanish, were 
they to be evaluated on Ai £ , by virtue of Lemma l3~Tl Since they are evaluated 
on JC £ q+1 \ instead, which is C(e 9+1 )-close to M. £ , these terms are 0(e g+2 ) and 
therefore of higher order for each of the blocks, recall Lemma 14.11 Thus, 

A ll) = R\tB sL [A f Q ( tlg] + R\^B^[A s Q ( flgl (4.24) 
A S = I^ ± B' ± [A f Q^ 1 g]+^Bf ± [A.Qi\gl (4.25) 
A (i = R^ ± B s± [A f Q^,g} + R 2 { f q) ± B^[A s Q^,g], (4.26) 

where the remainders of 0(e q+2 ) have been ommited for brevity. Recalling the 
definition of the Lie bracket, we rewrite Eq. (J4.24|) as 

A,',') = ^^((D 9 )A,Q^-^(A f Q^ 

+ R\£B^ ([Dg)A,Q'S - j t (A.Q<$) ) , (4.27) 

where we recall that all of the quantities are evaluated at (y,ipr q+ i),e). Next, 
(Dg)A s and the two time derivatives in Eq. ()4.27|) are zero to leading order 
by Lemma [A. II and [20| Lemma A. 2], respectively. Therefore, to leading order 
Eq. ()4.27|) becomes 

Ago) = R\^B s Q \Dg) A)Q^\ (4.28) 
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Here, AP^ stands for the leading-order term in the asymptotic expansion of 
^iq)(y> V'Cg+i)^/)) an d the right member is the leading order term in the 
asymptotic expansion of (R 1 ^B s± (Dg)AfQ^j)(y,h £ (y),s). 

We derive a similar formula for M qq y First, we rewrite Eq. ()4.25|) as 

A S = R^{{D 9 )A f Q^-^{A f Q^ 

+ R\^B^ ([Dg)A s Q% - ± (A S Q^ . (4.29) 

Next, Q§> = 0(e*), Q$ = 0(1), R^f- = (9(1), and R 1 ^- = 0(e«), by the 
induction hypothesis (|4.8|) - (|4.9|) . Thus, (201 Lemma A. 2] implies that the two 
terms in Eq. (|4.29|) involving time derivatives are 0(e q+1 ) and therefore of 
higher order. Also, (Dg)A s is zero to leading order by Lemma IA. 11 and thus 

A &r) = R\l%B^{Dg),A)Q^\ (4.30) 

We now substitute A^ 0) and h¥s from Eqs. (jOHj) and (OOj) in the 
expression U^ q ) = (A^ ^) _1 A^ (? - ) to find the desired expression for U^ q ) in 
terms of Q^ q \ 

u iq , g) = (Q[f ) y 1 Q^ ) . (4.31) 

We also need an expression for U^ q ) in terms of blocks of R( q ), which we 
will use in Eqs. ()4.19|) - (j4.22|) . Since R^ has the near block-diagonal structure 
given by the induction hypothesis (|4.8j) - (J4.9|) and is its inverse, we find 

to leading order for each of the blocks and for q — 1,2, ... . Equations (J4.31j) 
and (|4.32|) lead to the desired expression for Ui q>q ^ in terms of R( q ), 

^) = -<iKo)) _1 - ( 4 -33) 

Step 3. We derive the leading-order asymptotics of the matrix L( q y 

Recall that = A^JAlK)' 1 . Moreover, by LemmaEHJ is strictly 
0(1) and Af\ is strictly 0(e q+l ). Hence, L {q) = L {q ^ q+1) e q+1 + 0(e q+2 ), with 
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= A(gg+i)A(go)) expression for A^ n was derived in Eq. ()4.28|) . 

so here we focus on A 2 qq+1 y 

Equation (|4.26|) and the definition of the Lie bracket imply that 

A?i = ^((D 9 )A fQ ^-±(A fQ ^ 

+ R^B^ ({pg)A.Q® - ± (a.q£)) . (4.34) 

Next, Q$ = 0(1), = 0{e q+1 ), Bffi- = 0(e q+1 ), and R$- = 0(1), by 
the induction hypothesis. Also, the time derivatives are 0{e) by [213 Lemma 
A. 2], and thus the two terms in Eq. ()4.34j) that involve time derivatives are 
0(e q+2 ). Last, (Dg)A s = 0(e) by Lemma IA~T1 Thus, we find 

a L + d = 4 s i + D<(^)o4^/ 0) - ( 4 - 35 ) 

Equations ()4.28|) and ()4.35|) yield the desired formula for L^g+i) in terms of 
?(<?), 

L(q,q+1) = A(g,g+1) ( A (g,0)) = ^(g,q+l) (^M)) ■ (4.36) 



the blocks of R( g ), 



Next, we recast Eq. ()4.36|) in terms of blocks of Q^ q \ in order to use it 
in Eqs. (I4.15|) - (j4.18j) . The matrix Rr q \ is the inverse of and has the near 
block-diagonal form given in (|4.9|) . Thus, 

H{q) -{-e^(QTr i Q^ +1 \QTr 1 (^fV 1 / ( } 

to leading order for each block and for q = 1,2, ... . Equations (|4.36|) and (|4.37j) 
lead to the desired expression for L( qA+ i) in terms of the blocks of Q^ q \ 



- {Q^y 1 Q { ?« +1) ■ (4-38) 



Step 4. We substitute the results obtained in Step 2 and Step 3 into the 
formulas (|4~T5)) - (fP2|) derived in Step 1. 

Equations ()4.15|) and (j4.18J) . together with the induction hypothesis and 
the estimates U {q) = 0{e q ) and L {q) = 0{e q+1 ), imply that Q^ +1) and Q|f s +1) 
remain 0(1). This concludes the estimation of these blocks. 

Next, we show that Q ( 2 q f +1) = 0(e q+1 ). First, Q ( 2 q f +1) and Q%j are equal 
up to and including terms of 0(e q ~ 1 ), by Eq. ()4. 16|) and the estimate on U( q y 
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Thus, Qy X ' % ^ — for i — 0, 1, . . . , q — 1, by the induction hypothesis on Qg/ • 
It suffices to show that Q^f~ = 0. Equation ()4.16|1 implies that 

Q^=Q^f-Q^U {m) . (4.39) 

The right member of this equation is zero, by Eq. ()4.3H) . and the estimation 
of Q^j" 1 ^ is complete. 

Finally, we show that <5^ +1 ' ) = 0(e q+2 ) to complete the estimates on 
the blocks of Q( q+1 K First, and Q±] are equal up to and including 

terms of 0(e q ), by Eq. (j4.17J) and the order estimates on U( q ) and L( q y Thus, 
Qis +1,4 ' > = for i = 0, 1, . . . , g, by the induction hypothesis on Q^J . It suffices 
to show that Q^ rl ' q+1 ' = 0. Equation (j4.17|) implies that 

n (q+l,q+l) _ n {q,q+l) n {q,0) j ( A Af\\ 

Vis "Wis +V2s (4.4UJ 

where the right member of this equation is zero by Eq. ([4.38)1 . The estimation 
of is complete. 

The blocks of Ri q ) may be estimated in an entirely similar manner, using 
Eqs. (j4~T3 HP2|) . instead of Eqs. (j4~TH jl - (j4^|l . and Eqs. KSty and (jOHl) . 
instead of Eqs. ()4.31|) and (|4.38|) . The proof of Theorem 13 .21 is complete. 



5 The Michaelis— Menten— Henri Model 

In this section, we illustrate Theorem 13.21 bv applying the CSP method to the 
Michaelis-Menten-Henri (MMH) mechanism of enzyme kinetics [T5t IT6] . We 
consider the planar system of ODEs for a slow variable s and a fast variable c, 

s 1 = e(-s + (s + K- A)c), (5.1) 
c' = s-(s + k)c. (5.2) 

The parameters satisfy the inequalities < e <^ 1 and k > A > 0. Only 
nonnegative values of s and c are relevant. The system of Eqs. ()5.1|) - ()5.2|) is of 
the form ()2.2)) - (j2.3j) with m — 1, n — 1, y — s, z — c, g 1 — — s + (s + k — A)c, 
and g 2 = s — (s + k)c. 
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5.1 Slow Manifolds and Fast Fibers 



In the limit as e j 0, the dynamics of the MMH equations are confined to the 
reduced slow manifold 

M o = {(c,s):c=-?—,s>0}. (5.3) 

S ~\~ K 

The manifold Ado is asymptotically stable, so there exists a locally invariant 
slow manifold Ai £ for all sufficiently small e that is 0(e) close to Ado on an Y 
compact set. Moreover, A4 £ is the graph of a function h £ , 

M £ = {(c,s):c = h £ (s),s>0}, (5.4) 

and h £ admits an asymptotic expansion, h £ = h + eh\ + e 2 h 2 + • • • . The 
coefficients are found from the invariance equation, 

s - (s + K)h £ (s) = eh' £ (s)(-s + (s + k - X)h £ (s)). (5.5) 

The first few coefficients are 

, s , , . kAs , , . kAs(2kA — 3As — us — k 2 ) 

hois) = — — , Ms) = t — ■ — -r, ftaf* _ 



(5.6) 

In the limit as e I 0, each line of constant s is trivially invariant under 
Eqs. (I5.1|) - (|5.2j) . These are the (one-dimensional) fast fibers with base 
point p = (s, h (s j) G A4 . All points on JFq contract exponentially fast to p 
with rate constant —(s + k). The fast fiber J-^ perturbs to a curve that 
is 0(e) close to T§ in any compact neighborhood of A4 £ . The fast fibers J 7 ?, 
p G .M £ , form an invariant family. 



5.2 Asymptotic Expansions of the Fast Fibers 

To derive asymptotic information about the fast fibers, we look for general 
solutions of Eqs. (|5.1|) - (|5.2|) that are given by asymptotic expansions, 

s(t;e) = £V Si (t), c(t;e) = ^^(t), (5.7) 

i=0 i=0 

where the coefficients and c, are determined order by order. 

Consider the fast fiber with base point p = (s, h £ (s)), and let (s A , c A ) 
and (s B ,c B ) be two points on it; let As(t) = s B (t) — s A (t) and Ac(t) = 
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c B (t) — c A (t). The distance between any two points on the same fast fiber 
will contract exponentially fast towards zero at the 0(1) rate, as long as both 
points are chosen in a neighborhood of M. £ . We may write 

As(t;e) = £VA Si (t), Ac(t;e) = £VAcj(t), (5.8) 

i=0 i=0 

where As^t) = sf(t) - sf{t) and Ac.it) = cf{t) - cf(t). The condition on 
fast exponential decay of As(t) and Ac(t) translates into 

A Si (t) = 0(e- Cat ), Aci(t) = C(e- Cct ), t -> oo, (5.9) 

for some positive constants C s and C c . We let (s A ,c A ) and (s B ,c B ) be in- 
finitesimally close, since we are interested in vectors tangent to the fast fiber. 



5.2.1 0(1) Fast Fibers 



Substituting the expansions ()5.7j) into Eqs. (J5.1j) - ()5.2|) and equating 0(1) 
terms, we find 

s' Q = 0, (5.10) 

c 'o = s -(s + k)cq. (5.11) 

The equations can be integrated, 

So(t) = s (0) = s , (5.12) 
Co(t) = -^+(c Q (0)-^-)e-^ t . (5.13) 

S + K \ S + Kj 

Hence, 

As (t) = As (0), (5.14) 

Ac (t) = Ac (0)e-( So+ ^ + (9 so Co(t))A So (0) + 0((A S o(0)) 2 ). (5.15) 

The points A and B lie on the same fiber if and only if 

As (0)=0. (5.16) 

Thus, Eq. (I5.15J1 simplifies to 

Ac (t) = Ac o (0)e-( so+K ) i , (5.17) 

and Aco(t) decays exponentially towards zero, irrespective of the choice of 
Aco(O). Hence, Aco(O) is a free parameter. 

We conclude that, to 0(1), any vector ( ^ J with a constant (a ^ 0) is 



tangent to every fast fiber at the base point. 
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5.2.2 0(e) Fast Fibers 



At O(e), we obtain the equations 

si = -s + (so + « - A) Co, (5.18) 

ci = si - (s + k)ci - sic . (5.19) 

Using Eqs. ()5.12j) and (|5.13|) . we integrate Eq. (|5.18|) to obtain 

*i(f) = Sl (0) - + S ° + / :~ A f c (0) - (1 - e-<«+*). (5.20) 



S + K S + « V S + K 

Therefore, at 0(e), 

Asx(t) = Asi(0) + S ° + * - A Ac o (0)(l - e -( S0+K )*). (5.21) 

S + K 

For the two points to have the same phase asymptotically, it is necessary that 
linit^oo Asi(t) = 0. This condition is satisfied if and only if 

Asi (0) = - So + K ~ X Ac (0). (5.22) 

S + K 

Next, Ci(t) follows upon integration of Eq. ([5. 19)1 . 

c x {t) = Cl (0)e-( so+K )* 

(-1(0) + S -^± (c (0) - -^-)) (1 - e -C**) 
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(«0 + K) 

O>(0) - U(0) + 3 -^± U(0) + ) te-C»«" 

Sq + kJ \ S + K V. S + K/ 

_ g ° + — ~ 9 A L(0) - - J2_V (e -2( S 0+«)t _ e -( S0 +«) t) 

(so + k) 2 V so + ^y ; 

2(s + k) V s + 



KXS V f„-(8Q+K)t 



v s + re) 



+ (s + «)t-l). (5.23) 



We infer from this expression that lim^oo Aci(t) = 0, as long as Eqs. (|5.22j) 
and ()5.16|) hold. Hence, Aci(0) is a free parameter, just like Aco(0), and the 
only condition that arises at 0(e) is (|5.22|) on Asi(0). 

We conclude that any vector 

.(-<}- 

P 



a 



° N + e -{ 1 -^) a I, (5.24) 
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with a and (3 constant (a ^ 0), is tangent to every fast fiber at the base point 
up to and including terms of 0(e). Any such vector may be written as the 
product of a free parameter and a constant vector (fixed by so), 



1 A 



{a + e(3)\ V «»W \+0(e 2 ). (5.25) 



5.2.3 0(e 2 ) Fast Fibers 



At 0(e 2 ), we obtain the equation 

4 = si(c — 1) + (s + K — A)ci. 
Direct integration yields 



(5.26) 



s 2 (t) = s 2 (0) + 



A 



co(0) - 



k(sq + K — A) 



k(s + k — A)(s + k — 2A) + A 2 s 



(s + k) 



+ A( f° + ^ A) ( c (0) 



(s + 
c (0) 



Sl (0) 



s 



S + K, 



+ 1 



2(s + «0 3 
A 



So 



s + « 



Cl (0) 



S + K 

kXso 

(So + «) 4 



kA 



+ 



(s + «) 
kA 2 s 
2(s + k) 3 



Si{0) + S + K-X 



S + K 



2s r 



S + K 



t 2 + ft(t), 



(5.27) 



where the remainder involves the functions e (s ° +K )* ) te <S ° +K ) i , t 2 e ( s °+' c ) t ) 
and e ^ 2 ( s o+«)*. From this expression we find 

As 2 (t) = As 2 (0) + (d S0 s 2 {t)) As o(0) + (d C0 s 2 (t)) Ac (0) 

+ (d si s 2 (t)) Asi(0) + («9 Cl s 2 (t)) Aci(0) + 0(2) + 0( e - Ci ),(5.28) 

for some C > 0. Here, d co is an abbreviation for the partial derivative <9 co (o), 
and so on, and 0(2) denotes quadratic terms in the multivariable Taylor ex- 
pansion. First, we recall that As (0) = by Eq. (|5.16|) . Next, we calculate 
the partial derivatives in each of the three remaining terms, 



d C0 s 2 {t) 



Asx(0) k(s + k — A)(s + k — 2A) + A 2 s 



(S + K) 



v s + K) 
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+ A( ;° + - : 3 A) fc (o) - -^-) - KX{ ; o+K ~ x) t, ^ 

(s + k) 3 V s + kJ (s + k) 3 

d Sl S 2 {t) = ; rr Co(0) 



so + k) 2 \ s + kJ (s + k) 3 

-j^w^ (5 - 30) 

d Cl s 2 (t) = 1 — • (5.31) 

We substitute these expressions into Eq. ()5.28|) . recall Eq. ()5.22|) . and carry 
out the algebra to obtain 

As 2 (t) = As 2 (0) + fl ^_)A Cl (0) 



+ T^+tf V l{ ) (so + nf 

+ 0(2) + 0(e- ct ), C>0. (5.32) 



In the limit t — ► oo, Eq. (|5.32|) yields the condition 
As 2 (0) = - (l - ) Aci(0) 



Finally, Ac 2 (i) vanishes exponentially, as follows directly from the conditions 
()5.22|) and ()5.33|) . Thus, no further conditions besides ()5.33|) arise at 0{e 2 ). 

We conclude that any vector 



a 



\ / - M-^- Q 

+ e V s °+v 

/9 



+ £ 2 f - ^ - ^) P - J^f (si(0) + K(so ^ + ^2 Aso ) a j _ (r) _3 j, 

7 



with a, /3, and 7 constant (a / 0), is tangent to every fiber at the base point, 
up to and including terms of 0(e 2 ). 
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5.3 CSP Approximations of the Fast Fibers 



We choose the stoichiometric vectors as the basis vectors, so 

A<°> = (4°>, ,f) =(J l), B m =(5) = (? I)- (5,5) 

The CSP condition B^g = is satisfied if c = h (s), so the CSP manifold 
tCe coincides with Ai . With this choice of initial basis, we have 

A (0) = W » = ( e "f^ ) • (5.36) 



5.3.1 First Iteration 



At any point (s, c), we have 



(1) = 


(!. 


'(1) = 


("A 



S + K- A/ -1 \ (!) / I 



+ e— — — ^ , ^= _ £=1 ), (5.37) 



(i) 4(1) 

22 5 /1 12 / ' 



^=(4?, (5.38) 



In the first step, we evaluate 4 and BjL on /Q to obtain 

Hence, the CSP condition, 

1 , , k(— s + (s + k — X)c) 
Bfag = s-(s + k)c- ^— — 2 >-L = , (5.40) 

\S -\- K) 



is satisfied if 



s k\s 2 k 2 Xs(s + k - A) o x 



S + K (S + K) 4 (S + K) 



Equation (jo~4l| defines /C^, the CSPM of order one, which agrees with M. t 
up to and including terms of 0{e)\ recall Eq. ()5.6|) . 



Then, in the second step, the new fast basis vector, A± , and its comple- 
-B(i), in the dual basis are evaluated on K.^, 



4 X) = ( ? J -e f *(*+La) J f k as(,+ k -a) J +C(£ 3 ), (5.42) 
V 1 / \ ( S + K )3 y \ ( S+K )6 y 

B 2 W = (4?, -4?) (5-43) 
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Thus, we see that A\ is tangent to the fast fibers at their base points up to 
and including terms of 0(e) as Eq. (|5.24j) (with a = 1, (3 = — ^j^^ ) implies. 

As a result, Cs approximates TT £ also up to and including terms of 0(e). 

Remark 5.1. If, in this particular example, one evaluates A± on JC^ as op- 
posed to K,^ as we did above, then the approximation of TT e is also accurate 
up to and including terms of 0(e). 



5.3.2 Second Iteration 



The blocks of A(i) are 



A 11 

A (i) 



A 12 

A (l) 



A 21 

A (l) 



A 22 
A (l) 



[s + K - A) 



(s + k) + e- 

S + K 
2 (C-1)(S + K-X) 



+ e 
s 

S + K 

e 2 



(s + k) 3 
c- 1 



( c -l) + (c--f_) 

S + K 



A(c - 1) + [(s + « - A)c - s] , (5.44) 



c + e- 



[s + k) 



A(c- 1) - [(s + k-X)c-s] 



(5.45) 



(c- 1)(s + k- A)(s + «-2A) 



+ A[(a + k - A)c - s] + (s + k - A) 2 c 



S + K 



(5.46) 



A(c-l) + (s + «-A)( 



+ r 



with remainders of 0(e 2 



S + K 

j(c-l)(s + «-A) 



(s + «0 J 



S + K 

A(c - 1) - [(s + k - A)c - s] , (5.47) 



In the first step, we update A 2 and Bh, and evaluate the updated quan- 



tities on K!i\ to obtain 



.4 



A 



(2) 
12 

(2) 
22 



1 + e- 



,k\(2s — k)(s + k — A) 



K 



(s + k) 2 
+ e 



+ e 



(s + k) 6 
k\(k — 3s) 
(s + k) 5 



(5.45 



2 k 2 X(7s - 2k)(s + k - A) + k\ 2 s(s - 2k) 



(s + k)* 



B (2) 



4(2) ,(2) N 
^22 ' ^12 / ) 



(5.49) 
(5.50) 
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up to and including terms of 0(e 2 ). 
The CSP condition 



B\ 2) 9 



, >. re(— s + (s + K - A)c) 

8 - (s + K C - £^ ^ ^ 

(s + re) z 

2 / (3s — re)(— s + (s + re — A)c) 
+ e reA 



(s + re) 5 

(2s — re)(s + re — A)(s — (s + re)c) 



(s + re) e 



+ 0{e c 



is satisfied if 



0. 



+ e- 



(5.51) 



reAs 2 kXs(2kX — 3As — res — re^ 



s + re (s + re 



(s + re) 7 



0(£ 3 ).(5.52) 



Equation f)5.52j) defines , the CSPM of order two, which agrees with M. e 
up to and including terms of 0(e 2 ); recall Eq. (|5.6p . 

Then, in the second step, we update and B?^ to obtain 



A 



(2) 
11 



S + K — X 
S + K 



— e 



s + re 



(s + re - A) 2 c 



s + re 



[(s + re- A)(s + re-2A)(c- 1) 



+ A[(s + re - A)c - s] 



(5.53) 



A 21 



l + e 



(s + re- A)(c-1) 



s + re 
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+ (s + re- A) c- 



s + re 
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s + re - A) (s + re - 2A)(c - 1) 
2s + k" 



2c 



5, 



(2) 



4(2) 4(2) N 
^21 5 ^11 / > 



s + re 



with remainders of 0(e 3 ). Evaluating these expressions on K,f \ we obtain 



(5.54) 
(5.55) 



A 



A 



(2) 
11 

(2) 
21 



s + re — A 9 re(s + re — 2A)(s + re — A) + A 2 s , . 
h£ — , (5.56) 



l-e 

+ e 



s + re (s + re) 4 

re(s + re — A) 
(s + re) 3 

2 (s + re - A)(re 2 (s + re - 2A) + reAs) + reA 2 s 



B {2) 



A m _ A {2)\ 



:s + re) f 



(5.57) 
(5.58) 
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with remainders of 0(e 3 ). Therefore, A\ is tangent to the fast fibers at 
their base points up to and including terms of 0(e 2 ), according to Eq. (|5.34|) 

(with a = 1,(3 = 7 = ('+^)^('+^y +lt ^ +wAa ' ), and Cf ] is an 

C(£ 2 )-accurate approximation to TT £ . 

Remark 5.2. If one evaluates, in this particular example, A± on K e instead 
of on K,^ as we did above, then the approximation of TT £ is also accurate up 
to and including terms of 0(e 2 ). 

6 Linear Projection of Initial Conditions 

The main result of this article, Theorem 1.1 2\ states that after q iterations 
the CSP method successfully identifies TT £ up to and including terms of 
0(e q+1 ), where this approximation is given explicitly by A[ . This information 
is postprocessed to project the initial conditions on the CSPM of order q. In 
this section, we discuss the accuracy and limitations of this linear projection. 

Geometrically, one knows from Fenichel's theory that any given initial 
condition xq sufficiently close to M. £ lies on a (generally nonlinear) fiber 
with base point p on M. £ . Hence, the ideal projection would be itf(xo) = p 
(the subscript F stands for fiber or Fenichel) and this is, in general, a nonlinear 
projection. 

Within the framework of an algorithm that yields only linearized infor- 
mation about the fast fibers, one must ask how best to approximate this ideal. 
A consistent approach is to identify a point on the slow manifold such that 
the approximate linearized fiber through it also goes through the given initial 
condition. This approach was used, for example, by Roberts ^7] for systems 
with asymptotically stable center manifolds, where we note that a different 
method is first used to approximate the center manifold. Also, this approach 
is exact in the special case that the perturbed fast fibers are hyperplanes which 
need not be vertical. In general, if Xq lies on the linearized fiber and if 
ttf(^o) — P2, then the error \\pi — jp 2 1| made by projecting linearly is 0(s) and 
proportional to the curvature of the fiber (see also [T7j). 

For fast-slow systems, there is yet another way to linearly project initial 
conditions on the slow manifold. One projects along the approximate CSPF 
to the space T V T £ , where p is the point on the CSPM that lies on the same 
e = fiber as the initial condition. This type of projection is also consistent, 
in the sense that it yields an exact result for e = 0, but has an error of 0{e) 
for e > 0. Moreover, it is algorithmically simpler, since it does not involve a 
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search for the base point of the linearized fiber on which the initial conditions 
lie. However, it has the disadvantage that the projection is not exact in the 
special case that the fast fibers are (non-vertical) hyperplanes. 



A The Action of the 0(1) Jacobian on T p A4q 

The spaces T v T e and T p M. e depend, in general, on both the point p G M. e and 
s. As a result, the basis A also depends on p and e, and hence Aj and A s 
possess formal asymptotic expansions in terms of e, 

i=0 1=0 

Next, we compute the action of the Jacobian on A s to leading order. 



Lemma A.l Ker[Dg(p))o = T p A4q, forp G A^o- In particular, (Dg)oA° s = 0. 



Proof. The Jacobian is a linear operator, so it suffices to show that every col- 
umn vector of a basis for T p Ai vanishes under the left action of the Jacobian. 

We choose this basis to be the matrix ( J. 171 , 

V D y h o 

We compute 

Dgo ( Dyho ) = ( D y g 2 D z g 2 ) ( D™h ) = ( D y g 2 + D z g 2 D y h ) ' (L2) 

Differentiating both members of the 0(1) invariance equation g 2 (y, ho(y), 0) = 
with respect to y, we obtain 

D y g 2 {y, h {y), 0) + D x g 2 (y, h {y), 0)D y h (y) = 0. (1.3) 

Equations (jl.2|) and f)l .3|) yield the desired result 

d *U:0 = (o)' wM °- (L4) 



Finally, the identity (Dg) Q A Q s = follows from the fact that A° spans 
T P M , since A° s = A s \ £=0 by Eq. JEH). I 
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